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1. Introduction 

This final technical report covers a three and one-half year period preceding February 28, 
1993 during which support was provided under NASA Grant NAG- 1-1065. Following a gen- 
eral description of the system identification problem and a brief survey of methods to attack 
it, the basic ideas behind the approach taken in this research effort are presented. The results 
obtained are described with reference to the published work (references [1-23] in Section 5), 
including the five semiannual progress reports previously submitted ([2,6,10,12,16]) and two 
interim technical reports ([7,17]). 


2. List of Scientific Collaborators 

Y. A. Fiagbedzi 
A. V. Fullerton* 

J. Q. Pan 
A. A. Pandiscio 
A. E. Pearson* 

Y. Shen* 


Assoc. Prof, of Math. Sci., King Fahd Univ. 
Grad. Student Research Ass’t. 1989-91 
Grad. Student Research Ass’t. 1989-92 
Raytheon Grad. Student Fellow 1990-93 
Professor and Principal Investigator 
Grad. Student Research Ass’t. 1991-93 


* Received partial support under NAG-1-1065. 

3. The System Identification Problem 

The identification problem for a system with an input u and an output y is easy to state 
in generic terms: Given a class of systems with parametrized members Sq such that an 
input/output pair (u,y) is modeled in the equation-error format S q(u ,y )=0, find 8 such that a 
suitably chosen norm measure of distance, ||Sg(u,y)||, is “small” in some sense relative to all 
available i/o pairs (u,y). The models considered in this research are those that can be 
described by input/output ordinary differential equations. In the general nonlinear system 
case, the parametrized members are represented implicitly by differential operator equations of 
the form: 

Sq : ^ ^gj(Q)F jk (.u,y)P jk (p)E k (u,y) = 0, g Q =l (1) 

j = o *=i 

where the gj (0) are specified functions of the parameters 0, the Pj k (p ) are specified polyno- 
mials of degree n in the differential operator p=d/dt, and the (E k (u ,y),F jk (u,y)) are 
specified functions of the i/o pair ( u,y ). Here the problem is to obtain an estimate of the 
parameters 0 given sampled versions of the i/o data [u (/ ),y (t )] over some time interval 
0 <t<T . Since parametrized state vector equations of the normal form: 

x — f (x ,M ,0), y = h(x,u ,0) (2) 

constitute the most common starting point for methods aimed at the parameter identification of 
nonlinear differential systems, the first step in our approach is to obtain (assuming this is pos- 
sible) an equivalent i/o model like (1) when starting from a model of form (2). Methods for 
attacking the parameter identification of the model (2) include: quasilinearization, nonlinear 
filtering or smoothing, extended Kalman filtering, and state variable filters. The results of our 
efforts pertaining to the model (1) will be summarized in Section 4.8. 
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The more specialized but nevertheless omnipresent linear system counterparts to (1) and 
(2) are, respectively: the i/o differential operator models of the form 1 

Lq : D (p)y (t) = N (p)u(t) (3) 

in which 0 is comprised of the coefficients entering into the definitions of the polynomial 
matrix pair (D(p)fJ(p)), and the linear state vector equation model 

x(t) = Ax(t)+Bu(t), y(t) = Cx(t). ( 4 ) 

In contrast to (1) and (2), the relations between (3) and (4) are well known, although the 
multi-input/multi-output case entails cumbersome details. Thus, if G(s) denotes the transfer 
function matrix: 

G(s) =D~\s)N(s) = C (si -A )~ l B. (5) 

System identification for the class L e will also focus on determining the frequency transfer 
function G (i co) given an i/o pair under transient conditions (Section 4.9). 

Notwithstanding the problem of estimating the parameter vector 0, there is also the prob- 
lem of order or structure determination for whatever class of models is under consideration. 
Effort in this regard has been confined to the linear class L e . In addition to the system 
identification approaches mentioned above for the general nonlinear model (2), the methods 
applicable to the models (3) and (4) include: approximation techniques like the 5 operator, 
block pulse and orthogonal function expansion methods; the various moment functional tech- 
niques; and the conversion of an identified discrete-time model (obtained by a variety of 
methods) into a continuous-time model using, for example, the bilinear transformation. 2 3 

4. Research Results Under NAG-1-1065 

The system identification research carried out under this grant has utilized the method of 
moment functionals which employs integration-by-parts to transfer derivatives on data vari- 
ables to smooth user-chosen functions in order to obtain an algebraic equation useful for 
identification purposes. The particular method examined is originally due to Shinbrot in 
which the user-chosen “modulating” functions possess the property that the resulting moment 
functionals are devoid of unknown initial/boundary conditions on the data, i.e., no state or ini- 
tial condition estimation is involved for transient time-limited data. As introduced by Shin- 
brot, 4»(r ) is a modulating function of order n on a fixed time interval 0 </<T if it is 
sufficiently smooth and satisfies the In end point conditions: 


1 Here we allow for multi-inputAnulti-output systems. 

2 Sources for discussion of these methods can be found in the recent books: System Modeling and 
Identification, by R. Johansson, Prentice Hall, 1993. Identification of Continuous-Time Systems, N. K. 
Sinha and G. P. Rao, eds., Kluwer Acad. Pub. 1991. Identification of Continuous Systems, by H. Un- 
behauen and G. P. Rao, North-Holland, 1987. 

3 Shinbrot, M., “On the analysis of linear and nonlinear systems.” Trans, of the ASME, 79 , pp. 
547-552, 1957. Also, NACA Rep. No. TN 3288, “On the analysis of linear and nonlinear systems 

from transient response data” by M. Shinbrot, 1954. 
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p k $(t) = 0 at t = 0 and t - T, k = 0,1 • • n-1 (6) 

where p denotes the differential operator, i.e., p=d/dt , p 2 =d 2 /dt 2 , etc. Moreover, 
m=l,2 • • M, is a sequence of order n modulating functions if (6) is satisfied for each m. 

Although Shinbrot’s method and its unique property of avoiding state estimation has 
been known for a long time, it has remained relatively obscure due largely to a lack of clear 
guidance as to a good choice of modulating functions. In addition, the method may lack 
appeal on the basis of it being “off-line”. As we reviewed in 1985, 4 the ‘cons’ have 
apparently outweighed the ‘pros’ since very few papers have appeared on the method subse- 
quent to Shinbrot’s original work. 5 What is needed is a sequence of modulating functions with 
highly composite members that ease the computational burden in calculating the underlying 
“moment functionals”. The trigonometric Fourier modulating functions introduced in Pear- 
son and Lee (Footnote 4) satisfied this need because the required funct ional s are Fourier series 
coefficients of the i/o data which can be calculated efficiently by DFT/FFT techniques. How- 
ever, they were not presented in an explicit form, i.e., an a priori finite set of such modulating 
functions could be calculated by inverting a Vandermonde matrix depending on the order n 
and the number M of functions desired, but increasing n and/or M required a new inversion. 6 7 
This implicit stigma was removed by a discovery (see Eqn. (7) below) made during the early 
portion of the research carried out under NAG- 1-1065. The explicit Fourier modulating func- 
tions facilitated several extensions of the theory and helped pave the way to more sophisti- 
cated methods for handling noise corrupted data. 

4.1. The Explicitly Defined Fourier Modulating Functions 
Consider the modulating function of order n defined by 

♦*.(') = 0<t£T = 2nAi> 0 (7a) 

- < 7b > 

1 j = 0 v 

where i=V^ T, co 0 =27t JT plays the role of a ‘resolving’ frequency, and m is any integer which 
shall be referred to as the ‘modulating frequency index’. The fact that § m (t) is a modulating 
function of order n for each m is evident from the basic definition (6) applied to (7a). Thus, 


4 A. E. Pearson and F. C. Lee, “Parameter identification of linear differential systems via Fourier 
based modulating functions,” Control-Theory and Adv. Tech., I , pp. 239-266, 1985. 

5 This remains true today. In fact, the recent book edited by Sinha and Rao (see Footnote 2), 
although a good source for a variety of topics dealing with differential system identification, contains 
no reference to Shinbrot’s method. 

6 In the paper: “System identification using modulating fiinctions and fast Fourier transforms,” by 
T. B. Co and B. E. Ydstie, Computer Chem. Engr., 14 . pp. 1051-1056, 1990, independent verification 
has been given in the use of the real Fourier modulating functions by application to a variety of ODE 
models, as well as showing how the method can be used in model reduction. 

7 Equation (7b) follows from (7a) upon utilizing the binomial expansion. 
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the (|> m in (7), m=0,l, • • M, represent an explicitly defined linearly independent set of order 
n modulating functions with the important property that moment functionals involve Fourier 
series coefficients on [0,7]. 

Two modulation properties can be developed that ease the various least squares formula- 

^ /I —k 

tions. The first concerns a linear differential operator Pip), i.e., a polynomial like X a *P 

k = 0 

and a function v(r) on [0,7] with Fourier series coefficients V{m) defined by 

V(m) = ^-jvit)e~ ima)of dt. (8) 

1 o 

Property 1. For a sufficiently smooth function v (r ) defined on [0,7] and a differential opera- 
tor P (p) of order n, the integration of § m (t)P ip)v(t) over [0,7] satisfies 

T 

J<p m (f )P ip )v it )dt = A n P ( im co 0 ) V (m ) (9) 

where A" is the n ,h order finite difference operator, i.e., with Q (m )=P iim co 0 )V (m ), 

&Qim) = Qim+l)-Qim) 

A 2 Q im ) = Q im +2) - 2Q (m +1 ) + Q im ) 

m 

( 

A "(2(171) = t(- 1)* b)Qin+m-k). (10) 

k=0 

Proof. Taking into account (6) and using integration-by-parts on the left side of (9), all boun- 
dary point evaluations are zeroed while transferring the derivatives on v (r ) to derivatives on 
the modulating function <]> m (0 such that: 

t t 

fd> m it )P ip )v it )dt = \v it )P i-p )<|> m it )dt. 

0 0 

In view of the representation (7b), coupled with the fact that Pi~p)e~ }a =P(k)e~ )a , the right 
sides of (9) and (10) are obtained which proves the property. 

The implication of (9) is that derivatives of signals in the continuous time domain can be 
traded exactly for finite differences of modulated Fourier series coefficients in the discrete fre- 
quency domain by forming moment functionals using the Fourier modulating functions. The 
presence of such Fourier series coefficients should not be interpreted as implying a kind of 
Fourier series approximation in the identification problem. Rather, Property 1 facilitates an 
equivalent statement of the system identification problem entirely in the frequency domain 
while utilizing transient i/o data. In applications, the range of values to be chosen for m and 
the choice of the time interval [0,7], equivalently the resolving frequency co 0 , will depend on 
the bandwidth of the signals involved, as well as the complexity of the model relating the sig- 
nals. It is noted that A n X can be calculated by the M-file ‘diff(X,n)’ in MATLAB, which 
facilitates ease of coding in that environment, and the Fourier series coefficients of the signals 
required in the identification problem can be calculated efficiently and accurately via 
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DFT/FFT techniques. Further elaboration on this point follows below in Section 4.3. 

4.2. Least Squares Identification for L e : SISO Case 

Let an i/o data pair [u(t),y(t)] for a single-input/single-output linear system of order n 
be modeled on [0,r] by 

p n y(t) + '£a j p n -jy(t) = Zbjp n - j u(t) + e(t) (11) 

y=i J= l 

where the (fl ; ,&/),;'= 1,2 ■ ■ n, are parameters and e(t) represents the effect of modeling 
errors. Modulating (11) with <t> m (f), integrating over [0 ,T] and applying Property 1: 

A" (im coq)" Y(m) + X^A" (im(d 0 ) n ~ j Y(m) = ^bjA n (imGi 0 ) n ~ J U(m) + e(m) (12) 

7=1 f=i 

where the residuals in (12) are defined by 

T 

z(m) = \e(t)$ m (t)dt. (13) 

o 

Defining a column vector 6=col[-<3 j, • • -a n ,b h • • b n ], a rearrangement of terms in (12) 
leads to the standard (albeit complex-valued) linear regression: 

yj(m) = 7(m)0 + e(w) ) m = 0, 1, • • M (14) 

where the row vector of regressors y(m ) is defined by 

7 (m) = row( y{(m), ■ ■ y%(m\y“(m), • • Y“( m ) ) ( 15 ) 

and the pairs (yj(m),yj(m)), y'=0,l ■ • n , are defined by 

[y f(m ), yj(m ) j = A" (im U> 0 ) n ~ j \u (m ), Y (m ) j . (16) 

Taking into account the fact that (14) is complex-valued, a cost function for least squares 
minimization can be defined by 

J(Q) = (Y C -r c ey(Y c -T C Q) (17) 

where prime denotes transpose and the following notation applies: 



Re yJfcO) 


Re Y(0) ' 
* 

Y c = 

ReY $(M) 
Im Y^(0) 

. Tc = 

• 

ReY (M) 

Im y(0) 
• 


ImYj(M) 


IrntfA/). 


with the complex y(m) and y$(m) defined in (15) and (16). 

Assuming linearly independent regressors, the one-shot estimate which minimizes (17) is 
given by 
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e = (r c T c r 1 iyy c . ( 19 > 

Noting from (10) that M in (18) actually corresponds to the frequency (M +n )co 0 suggests that 
in choosing the pair {M ,C0q), their product should correspond roughly to the system 
bandwidth. Also, M has to be sufficiently large to provide enough algebraic equations upon 
which to base a one-shot estimate. Hence, if # denotes the number of unknowns (# = 2 n for 
the model (11)), then the pair (M,co 0 ) can be uniquely determined by the guidelines: 

M ~ 2#~ 4# 


and ( M+n )CDo - ( 0 # (20) 

where tOg is the system bandwidth, i.e., the highest useful frequency in the i/o data. 


4.3. Computational and Sampled-Data Considerations 


Assuming uniform sampling of the i/o data on [0,7’], the integral (8) for all required 
Fourier series coefficients will be approximated by standard quadrature formulae which, in 
turn, are calculated by a DFT/FFT algorithm. Thus, if the sampling interval 8 1 and the total 
number of samples A satisfy 5r=77A and N>M+n , the staircase, trapezoidal and parabolic 


rules applied to (8) lead respectively to 


■ jj 


A 
v 0 +v N 


\N - 1 
k=0 


mk 


+ 0 <N> 


N - 1 

+ 

*=i 


mk 


1 


+ 0 H 7 ) 

N 2 


(21a) 

(21b) 



-im OV 



VQ+Vy 

2 


N - 1 . N - 1 . 

+ L 2v k W mk + X ^>k w 
it =1,3 • • *=2,4 ■ • 


+ 0(-t) (21c) 
A 4 


where 

m = 0, 1, • • M+n, v* =v(kT/N) and W = 

The quantities within brackets on the right hand sides of the above can be evaluated by a 
standard DFT/FFT algorithm provided the range of values for the modulating frequency index 
m is artificially extended to [0, A-l]. A fundamental property of the DFT is that the first 
and last halves of the transformed sequence are complex conjugates of one another. Hence, in 
order to preserve uniqueness for the first M+n harmonics, it follows that N has to satisfy 

A/2 > M + n. (22) 


In practical experience using either simulated or physical data, very little difference is 
noticed between the three approximations in (21) provided: (i) the data is oversampled in 
terms of the Nyquist sampling theorem, and (ii) the guideline (20) is upheld. 8 These provisos 

8 The first proviso is more essential than the second if the “adaptive” weighted least squares 
(AWLS) algorithm is utilized (Section 4.4), i.e., the AWLS algorithm tends to make the parameter esti- 
mates insensitive to a mismatch in bandwidth. 
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also bear a relation to the DFT requirement (22) as can be shown by the following. Let F B 
and F s denote the bandwidth and sampling frequency respectively (in Hertz), then ‘oversam- 
pling’ implies the inequality: 


Assuming adherence 


2 F b <<F s 


K 


_1 

St 


<< — = — = 


N 

T 


A/co c 

2k 


of the bandwidth guideline (20), the above implies 


(Af +n )Oq N (0 0 


< < 


71 


2n 


M_ + n 

N/2 


<< 1. 


(23) 


Hence, the DFT requirement (22) is amply upheld by oversampling and adherence to the 
bandwidth guideline in (20); moreover, (23) implies that only the lowest few percent of the 
available DFT harmonics are utilized in the least squares identification. 

The following example from [22] illustrates some aspects of the sampling rate considera- 
tion and, at the same time, provides a comparison with a commercially available algorithm for 
parameter identification. Consider the system 

y(t) + 3 y(t) + 8 y(t) = 5*1(0, 0 <t<T = 10s (24) 

which possesses the transfer function: H(s) = 5/(s 2 +3s+8). Shown in Fig. 1 for point of 
reference is: (a) the time domain step response of the system, (b) the frequency domain mag- 
nitude plot, and (c) the effect of sampling rate on the parameter estimation errors under ideal 
zero-noise conditions. The latter is discussed below. Since there are only 3 unknown param- 
eters and n =2, M was chosen as Af= 6. This means that (A/+n )CDo=0.8 Hz is the highest fre- 
quency that will be extracted from the i/o data, which is roughly the bandwidth of the system 
as seen in Fig. 1(b). 
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Fig. 1. Aspects of the simulated system 
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The output y(t ) was corrupted by additive white Gaussian measurement noise. Two 
hundred Monte Carlo runs were made for each of several noise-to-signal ratios (NSR) under 
two separate conditions: (i) the initial conditions fixed at (0,0), and (ii) the initial conditions 
randomized for each run, i.e., a total of 400 Monte Carlo runs at each NSR. The input was 
K(f)=sinf 2 /5 over the T =105 time interval for each run. All calculations were earned out in 
MATLAB. The sampling rate was fixed at 25.6 Hz, thus facilitating a standard 256 point 
DFT for the 10s time interval, and the staircase rule (21a) was used in approximating (8) to 
obtain the i/o Fourier series coefficient pairs [U (m),T(m)], m=0,l • • M+n. Since N = 256 is 
much larger than A/+n=8, it is noted from (23) that only 8/128 ~ 6% of the total harmonics 
available from the i/o DFT’s is actually retained in implementing the algorithm. This 
accounts for the high accuracy achieved in implementing the method with modest sampling 
rates inasmuch as calculating the higher frequency coefficients would be expected to entail 
greater numerical error in the presence of high frequency noise. This is further discussed 
below. 






(a) NORMALIZED BIAS (%) 
VERSES NSR (%) 


(b) NORMALIZED STD (%) 
VERSES NSR (%) 


Fig. 2. Measurement noise effects 
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The results are shown in Fig. 2 where “MFT” means the estimate is based on the modu- 
lating function technique. Also shown is the estimate based on the prediction error method 
(PEM) from the Identification Toolbox in MATLAB. 9 In the case of nonrandom (fixed) initial 
conditions shown in the top halves of Fig. 2, the PEM shows a fairly constant 6=7% normal- 
ized bias error 10 as the NSR increases from 0% to 60%, whereas the bias for the MFT gradu- 
ally increases from 0% to =12% over the same range of NSR values. Meanwhile, the normal- 
ized standard deviations are nearly the same for this case. However, in the case of random- 
ized initial conditions, see the lower halves of Fig. 2, both the bias and the standard deviation 
for the PEM show a several-fold increase and, moreover, each exhibits rather unpredictable 
behavior that is difficult to reproduce even with additional Monte Carlo runs. By contrast, the 
bias and standard deviation have both decreased at each NSR for the MFT, indicating that 
nonzero initial conditions favor the MFT. 

The above comparative results have been verified on other examples as well (see Fuller- 
ton [11]). One reason for this is that the MFT does not have to estimate unknown initial con- 
ditions; another is related to the fact that the MFT is a direct identification technique for 
continuous-time models, while the PEM first estimates the parameters of a discrete-time 
model then converts this to a continuous-time model. Also, the Fourier series coefficients 
needed by the MFT can be computed accurately with modest sampling rates using the DFT. 
Additional insight in this regard is gained by referring to Fig. 1(c) which shows the influence 
of sampling rate on the percent normalized error under zero measurement noise conditions. 
Apparently the PEM needs a much higher sampling rate for the given 10s of data, or a much 
longer T time interval for the given 25.6 Hz sampling rate, in order to achieve a small estima- 
tion error. The influence of sampling rate on the estimation error for the MFT is much more 
consistent with the Nyquist sampling theorem in view of the frequency magnitude plot in Fig. 
1(b). Finally, it is interesting to note that the MFT requires substantially less computer time 
for each set of Monte Carlo runs (about a ten to five fold decrease depending on whether or 
not the PEM has to estimate the initial conditions). 

4.4. The LS, WLS and AWLS Algorithms 

If the equation error function e(t) in the linear system model (11) is modeled as station- 
ary Gaussian “white” noise with zero mean and covariance Ee (t)e(t+x)=ob( t), then the 
modulated residual process e(m) defined in (13) is also Gaussian with a covariance matrix 
given by Equation (8) in [16], or (2.25) in [20]. Even though the residuals in the 
identification problem will not be white, this fact facilitates a reasonable weighting matrix in 
the discrete frequency domain for a weighted least squares estimate which attempts to 
ameliorate the effects of noise when the measurement noises on the i/o data are white. More- 
over, it points the way to formulating an “adaptive” weighted least squares algorithm in the 


Ljung, L. (1987). System Identification Toolbox, Version 2.11. The MathWorks, Inc. 


10 Defined as 
mates for the true Qj * . 


1 # 
jx 

ff y=i 


e,-V 


e ; * 


• 100% where 6j is the ensemble average of the MFT esti- 
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discrete frequency domain which emulates an approximate maximum likelihood estimate. 
These are the basic ideas developed in Chapter 2 of Shen’s Ph.D. thesis [20]. The advantage 
of these developments is illustrated by a simulation study carried out for the system (24) and 
summarized in Section 2.4 of [20]. In addition to the sampled i/o data on [0,T ] , implementa- 
tion of the AWLS algorithm requires an assumption of the system order n and an initial guess 
for the weighting matrix W which is normally taken as the covariance matrix under the ideal 
white noise residual case. Convergence is fast and has never failed in numerous trials. The 
software for these algorithms will be supplied under separate cover in a toolbox to be 
prepared by Y. Shen and the PI. An earlier toolbox based on the real Fourier modulating 
functions was prepared by Fullerton [8]. 

4.5. Model Reduction 

Motivated in part by the requirements attending one of the two methods developed in 
this research to handle MIMO models, the AWLS algorithm was applied to the model reduc- 
tion problem, thus giving a frequency domain approach to this problem that has several 
advantages: (i) the high order model to be simplified can be supplied in either the time 
domain via a differential equation, or in the frequency domain via Bode type plots, (ii) 
weighting the cost function in the frequency domain is a straightforward and direct option for 
the AWLS formulation, and (iii) the resulting lower order model obtained via the AWLS for- 
mulation is at least as good and usually better than the same order model obtained by all 
methods with which the algorithm has been compared. Most notably, this is true for the 
well-known Balanced Reduction method, which is available in MATLAB. The detailed dis- 
cussion and verification of this is given in Chapter 3 of Shen’s thesis [20]. 

4.6. Extensions for L 0 : MIMO Systems 

Two viewpoints in modeling multi-input/multi-output systems led to two different formu- 
lations to extend the LS, WLS and AWLS algorithms for the class L 0 in (3). If the underly- 
ing model is selected to be a state space model as in (4), then an appropriate scalar-valued i/o 
differential operator model, which is required in this approach, is the following: 

m u 

P(p)yj(t)= ZQjk(P)u k (t), 0 <t<T (25) 

Jfc=l 

j = 1,2 • • m y 

where the scalar-valued polynomial P(s ) is defined by P (s )=det( 5 '/ —A ) in relation to the sys- 
tem matrix A in (4). A joint cost function that reflects the equation errors for each equation 
in (25) can be specified and minimized using the AWLS algorithm. A weighting between 
each such equation error has to be supplied, or determined by testing different weightings, 
which essentially accounts for differences in the output measurement noise variances for the 
my outputs. This is the viewpoint taken in Chapter 5 of Shen’s thesis [20] (also summarized 
in [16]) which was used to model the physical data supplied under NAG-1-1065 for an F-18 
aircraft. 

Another viewpoint in modeling MIMO systems of the class L 0 is to find the best linear 
transfer function model that relates the i/o data without constraining the polynomials on the 
left hand side of (25) to be the same for each j. That is, each output y ; could be modeled 
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by: 

m u 

P,(p)y.U) = I,Q lt <P)u t <‘>, 0£1<7- (26) 

* = 1 

which differs from (25) in that the poles of the estimated MISO models are not constrained to 
be the same. This leads to a simpler extension of the AWLS algorithm developed for the 
SISO case, but necessitates a model reduction for each resulting MISO model in (26) in order 
to ferret out the common poles and zeros. This viewpoint is developed in Chapter 4 of 
Shen’s thesis [20] and motivated the model reduction results mentioned in Section 4.2 above. 

4.7. Applications to an F-18 Aircraft 

Using physical data for an F-18 aircraft supplied under NAG-1-1065, the AWLS algo- 
rithm was applied to obtain i/o models for both the lateral and longitudinal dynamics 
[16,17,20]. The following output signal-to-error ratio (in decibels) was used as a measure of 
goodness-of-fit in each channel: 

S/E = 20 logio ^i , eit) =y(f) -y(r), 0 <t<T 

where RMS(y) means the root-mean-square value of the output signal y(t) over the [0,7] 
time interval, and y {t ) is the estimated output using the model. Since our algorithm does not 
estimate initial conditions, we used a standard Luenberger observer running backwards in time 
(see details in Appendix C of Shen [20]) in order to estimate an initial condition and, subse- 
quently, simulate the model with the given input data to obtain e(t ) for the above S/E. Typi- 
cal S/E ratios were on the order of 8 - 12 db for the longitudinal dynamics, depending on 
whether a second or third order model is used, and 9 — 12 db for the lateral dynamics. (See 
Figs. 5.5 and 5.12 in Shen [20] for the time domain performance of the resulting models.) In 
order to compare with a simulated model under specified measurement noise conditions, the 
same algorithm was applied to a theoretical linear model supplied to us under NAG- 1-1065 
using the given forcing function data for the physical modeling. These models gave S/E 
ratios on the order of 34 ~ 41 db for the longitudinal dynamics and 16 - 30 db for the lateral 
dynamic channels (see Figs. 5.8 and 5.14 in [20] for the time domain performance of these 
models). Even at 16 db the fit is very good, and at S/E ratios above 30 db there is hardly any 
discernible error in the model responses compared with the theoretical outputs. 

4.8. Extensions to LTV and Nonlinear System Models 

If the time variations in the coefficients of a linear time-varying differential equation 
model can be parametrized as linear combinations of known modal (smooth) functions on 
[0,7], then the least squares formulation developed in Section 4.2 can be extended to yield 
modulated equation errors in the discrete frequency domain with Fourier series coefficients 
computed for time-modulated i/o data, i.e., the i/o data modulated by the modal functions on 
[0,7]. This is carried out in [22] and is illustrated by simulations for a second order system 
with a parabolically varying damping coefficient. The good noise rejection properties experi- 
enced with the algorithm for time-invariant systems appear to carry over to the time-varying 
case as well. 
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If the form of the nonlinearities are known up to parametric representation and if the i/o 
model can be manipulated into the generic form (1), then it is possible to extend the least 
squares parameter identification in the discrete frequency domain. This is earned out in a 
series of publications in [1,18,19]. Although mainly of theoretical interest thus far, the poten- 
tial exists to apply these algorithms to unsteady aerodynamic models akin to those studied in 

Chin and Lan. 11 

4.9. Nonparametric Identification: Frequency Analysis 

Methods for determining the frequency transfer function of a stable linear system from 
i/o data include correlation and spectral analyses, as well as the direct ratio of Fourier 
transforms and steady state sinusoidal measurements. Each of these nonparametric 
identification techniques require either a statistical stationarity assumption on the data, or a 
periodic steady state condition to be established, before initiating calculations of the transfer 
function at pertinent frequencies. Notwithstanding noise considerations, long data lengths may 
be required in order to achieve good accuracy due to the stationarity or steady state assump- 
tion, or to eliminate end point effects in case a direct ratio of Fourier transforms is used on 
time-limited data. By contrast, a method was proposed in [13,14] that utilizes the modulating 
functions (7) to extract the frequency content in short data lengths in order to set up a least 
squares estimation of the transfer function at selected frequencies. Since short data lengths 
are used there is no assumption of steady state operation or stationarity of the data, though 
there must be present sufficient energy content in the data at the specified frequencies in order 
to avoid degeneracy in the least squares estimate. The method was also extended to mul- 
tivariable systems, and was illustrated by way of simulation for a two-input/two-output exam- 
ple in [22]. 

4.10. Miscellaneous Notes 

In addition to the results summarized above, research under NAG- 1-1065 has included 
the following topics. 

Order Determination for L e The problem of determining the orders of the various polynomi- 
als that enter into specifying a member of the linear class L 0 has been studied in [9,14,20]. 
Basically, the approach taken is to examine the residuals £{m ) in the discrete frequency 
domain for various presumed orders starting from a lowest order. Since there are no initial 
conditions to be estimated in the modulating function approach, the norm of the residuals in 
the frequency domain is unbiased by unknown initial/boundary conditions and provides a use- 
ful measure for order determinaion. However, more could probably be done on this important 
problem and, hence, this remains a topic for future investigation. 

Collinearity The problem of collinearity in the regressors brought about by linear feedback 
was studied in [11]. The result is that predicting the deleterious effects of noise and collinear- 
ity remains a difficult problem in system identification, but the degree of collinearity that can 


11 Chin, S. and C. E. Lan (1991). Fourier Functional Analysis for Unsteady Aerodynamic Model- 
ing. AIAA-91-2867-CP. 
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be tolerated in the modulating function approach is quite high, e.g., about 95% degree col- 
linearity (as defined in [11]) while maintaining about 5% accuracy in the estimated parameters 
and in the presence of about 5% noise on the data. 

HO Modeling to Alleviate Noisy Measurements A study was undertaken in [7] to derive an 
i/o model not involving a noisy velocity signal with the expectation that the accuracy in the 
parameter estimates might be enhanced if the signal in question was particularly noisy. The 
state equation model with appropriate noise levels was communicated to us from E. A. 
Morelli. 12 The expectation was realized for those parameters that could be estimated from the 
i/o model, i.e., not all the state equation parameters could be estimated uniquely from the i/o 
model due to an inherent lack of identifiability for certain state equation parameters in the 
input/output model. 

Frequency Estimation in a Signal Processing Problem The classic problem of estimating the 
frequencies in a sum of sinusoids buried in noise can be formulated as a system identification 
problem of the class L 0 using the fact that the sinusoids satisfy a homogeneous linear 
differential equation. The modulating function approach was used to study this problem with 
the expectation that the avoidance of initial condition estimation could enhance the accuracy 
in the estimates. The results of the study are contained in [14,15] and basically confirm this 
expectation in relation to the high order Yule-Walker approach. 

Time Lag Systems Although not funded by NAG- 1-1065, research results have been obtained 
on an approach to the feedback control of time lag systems that employs a ‘ ‘reducing transfor- 
mation” [3-5,21]. The main advantage of this approach is that it facilitates the design of 
feedback controllers for systems described by differential-delay equations based on well- 
known finite dimensional control techniques. Our recent studies in [21] extend robust control 
techniques, specifically the structured singular value, to this class of control systems within 
the context of the reducing transformation method. 
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